'''
WARNING:
  Preserve Python3 compatibility
  print()
'''


'''
Jorg Stelling model from
A synthetic tunable mammalian oscillator
'''

import sys
from numpy import *
from scipy import *
from bricks import *


'''
ARGUMENTS
  TIME : simulation time (minutes)
'''
TIME = array(range(int(sys.argv[1])))


'''
MODEL PARAMETERS


G1: [57][59][pCMV][TAL97:KRAB][NEPTUN]
G2: [59][97][pCMV][TAL57:KRAB][mCITRIN]
G3: [PIP][pCMV][TAL97:KRAB]
G4: [RIP][pmin][TAL57:KRAB]
G5: [pCMV][Pristinamicin IP]
G6: [pCMV][Rapalog IP]
'''

#Gene dosage plasmid [plasmids/cell]
#Relative number to cell 
#TODO: to use with measurements transfer function
#TODO: HEK 293
#Batard 2001

#*Assumption: uniform distribution
#*Tunable
G1 = 200.0
G2 = 200.0
G3 = 200.0
G4 = 200.0
G5 = 200.0
G6 = 200.0



#Experimental TAL Effector rates
#RLA Measurements
tal97 = Effector('TAL97:KRAB', 'dat/tal97krab.dat')
tal57 = Effector('TAL57:KRAB', 'dat/tal57krab.dat')
#tal97.plotTransferFunction()
#tal57.plotTransferFunction()
#show()
#quit()
r97 = (G3+G1)/(G2)
r57 = (G2+G4)/(G1)
q97 = tal97.getTransferRate(r97)
q57 = tal57.getTransferRate(r57)
#print q97
#print q57
#quit()


#Inducer affinity
#(testing)
#TODO: look for source
Kpc = 1.0
Krg = 1.0

#Maximal transcription rates (mRNA/min)
#apply pCMV, pMin
#TODO: Jan!
#Apply protein length
k1 = 30.0
k2 = 30.0
k3 = 30.0
k4 = 30.0
k5 = 30.0
k6 = 30.0


#Basal activity of inducible promoter
#TODO: source, RLA, mammalian
#TODO: Jan, Miha
alfapc = 0.085
alfarg = 0.056

#Binding affinity 
#Literature data. PIP repressor.
#Tigges.
KPC = 3.0
KRG = 3.0

#Cooperativity
#TODO: source
npc = 2.0
nrg = 2.0

#Translation rates (protein/min)
#TODO: source
#TODO: protein length adjust
t57 = 0.02
t97 = 0.02
tpc = 0.02
trg = 0.02


#Degradation rate (mRNA/min)
#TODO: source
kdm = 0.0173

#Degradation rate (protein/min)
#TODO: source
kdp = 0.0058

#Concentration scaling factor for copies/cell
#TODO: source
fv = 1.85e-3


#Temporary switch
#TODO: sources
r1 = k1*fv
r2 = k2*fv
r3 = k3*fv
r4 = k4*fv
r5 = k5*fv
r6 = k6*fv




#TODO: Inducible repressor asymmetry
#sources
ka = 1/15.0






#################
# SIMULATION PART
#################



def plotResults():
  if sys.version_info[0] > 2:
    return

  from pylab import plot
  from pylab import show
  from pylab import legend
  from pylab import figure
  from pylab import semilogy
  
  f1 = figure()
  plot(TIME, RBSm, linestyle='dashed', label='RBSm', color='blue')
  plot(TIME, RBSp, label='RBSp', color='blue')
  plot(TIME, PBSm, linestyle='dashed', label='PBSm', color='green')
  plot(TIME, PBSp, label='PBSp', color='green')
  legend(loc="upper right")

  f2 = figure()
  plot(TIME, PBSa, label='PBSa', color='blue')
  plot(TIME, RBSa, label='RBSa', color='red')
  plot(TIME, PC, label='Pristinamicin', color='blue', linestyle='dotted')  
  plot(TIME, RG, label='Rapalog', color='red', linestyle='dotted')
  legend(loc="upper right")

  f3 = figure()
  plot(TIME, T97m, linestyle='dashed', label='T97m', color='blue')
  plot(TIME, T97p, label='T97p', color='blue')
  plot(TIME, T57m, linestyle='dashed', label='T57m', color='green')
  plot(TIME, T57p, label='T57p', color='green')
  legend(loc="upper right")
  #semilogy()
  show()
  
    

'''
OBSERVED SPECIES
'''
RBSm = array([0.0 for i in TIME]) #rapalog binding protein mRNA
RBSp = array([0.0 for i in TIME]) #rapalog binding protein
RBSa = array([0.0 for i in TIME]) #rapalog binding protein (active)

PBSm = array([0.0 for i in TIME]) #pristinamicin binding protein mRNA
PBSp = array([0.0 for i in TIME]) #pristinamicin binding protein
PBSa = array([0.0 for i in TIME]) #pristinamicin binding protein (active)

PC = array([0.0 for i in TIME]) #pristinamicin
RG = array([0.0 for i in TIME]) #rapalog

RG[int(3*len(TIME)/8.0):int(7*len(TIME)/8.0)] = 1500
PC[int(3*len(TIME)/8.0):int(7*len(TIME)/8.0)] = 1500

T97m = array([0.0 for i in TIME])
T97p = array([0.0 for i in TIME])
T57m = array([0.0 for i in TIME])
T57p = array([0.0 for i in TIME])

'''
SIMULATION
'''

for t in TIME:
  #Skip first step
  if t == 0:
    continue

  
  #TODO: check with article
  dPBSm = G5 * r5 - kdm * PBSm[t-1]
  dPBSp = tpc * PBSm[t-1] - kdp * PBSp[t-1]
  # dPBSa = PBSp[t-1] * (1 - PC[t-1]/(Kpc + PC[t-1]))
  dPBSa = PBSp[t-1] * (PC[t-1]/(Kpc + PC[t-1]))

  
  #TODO: check with article
  dRBSm = G6 * r6 - kdm * RBSm[t-1]
  dRBSp = trg * RBSm[t-1] - kdp * RBSp[t-1]
  #dRBSa = RBSp[t-1] * (1 - RG[t-1]/(Kpc + RG[t-1]))
  dRBSa = RBSp[t-1] * (RG[t-1]/(Krg + RG[t-1]))

  r3 = fv * ka * k3 * (alfapc + (1-alfapc)*(PBSa[t-1]**nrg/(KPC + PBSa[t-1]**npc)))**-1
  r4 = fv * k4 * (alfarg + (1-alfarg)*(RBSa[t-1]**nrg/(KRG + RBSa[t-1]**nrg)))
  
  
  s = ''
  
  #TAL 97
  dT97m = G1 * r1 - q57 * T57p[t-1] 
  T97m[t] = T97m[t-1] + dT97m
  T97m[t] = int(T97m[t] > 0) * T97m[t]
  s = s + str(T97m[t]) + ' ' 
  dT97m = G3*r3 - kdm * T97m[t-1]
  T97m[t] = T97m[t] + dT97m
  dT97p = t97 * T97m[t-1] - kdp * T97p[t-1]

  
  #TAL 57
  dT57m = G2 * r2 - q97 * T97p[t-1]
  T57m[t] = T57m[t-1] + dT57m
  T57m[t] = int(T57m[t] > 0) * T57m[t]
  print s + str(T57m[t])
  dT57m = G4*r4 - kdm * T57m[t-1]
  T57m[t] = T57m[t] + dT57m
  dT57p = t57 * T57m[t-1] - kdp * T57p[t-1]


  #print t, r3, r4,  dT57m, dT57p, G4*r4


  PBSm[t] = PBSm[t-1] + dPBSm
  PBSp[t] = PBSp[t-1] + dPBSp
  RBSm[t] = RBSm[t-1] + dRBSm
  RBSp[t] = RBSp[t-1] + dRBSp   
  
  PBSa[t] =  PBSa[t] + dPBSa
  RBSa[t] =  RBSa[t] + dRBSa
  
  
  T97p[t] = T97p[t-1] + dT97p
  T57p[t] = T57p[t-1] + dT57p

  #Preserve above zero
  PBSm[t] = int(PBSm[t] > 0) * PBSm[t]
  PBSp[t] = int(PBSp[t] > 0) * PBSp[t]
  RBSm[t] = int(RBSm[t] > 0) * RBSm[t]
  RBSp[t] = int(RBSp[t] > 0) * RBSp[t]

  PBSa[t] = int(PBSa[t] > 0) * PBSa[t]
  RBSa[t] = int(RBSa[t] > 0) * RBSa[t]

  
  T97p[t] = int(T97p[t] > 0) * T97p[t]
  T57p[t] = int(T57p[t] > 0) * T57p[t]


#plotResults()
